% GAUTIER LE BIHAN - 2020
% Replication files for “Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model�? –Review of Economics and Statistics
%
% Compute additional simulated moments for Appendix figure H
% Investigate identification based on the Canova & Sala approach
% Relying on codes by Nakamura and Steinsson 
%**************************************************************

clear;
tic;
disp('Program started running on date/time:');
disp(datetime);

addpath('..\..\Utilities')  



% Set sector specific parameters
%***************************************************************
% NB: Edges of the rad grid need to be extra padded in the multi-sector 
% model because the sectors have different mean rad.

NSectors = 2;
sm = 0.0;

%0.0552582422061334 0.0345297530098889 -3.22066898186641 0.850810309615751


param0=[0.0521    0.0348   0.0486   0.60];
param0=[0.0503 0.0303 0.0344   0.6946];
%[0.0521115504741103 0.0348316791205970 0.0486324094042507 0.538186216692536 0.0166306069512913]
s=10000;
p0r = unifrnd(param0(1)-0.005,param0(1)+0.005, s, 1) ;
mur = unifrnd(param0(2)-0.01,param0(2)+0.01, s, 1); 
sigr = unifrnd(param0(3)-0.01,param0(3)+0.01, s, 1); 
rhor = unifrnd(param0(4)-0.1,param0(4)+0.1, s, 1); 
param=[param0; p0r mur sigr rhor];


s=s+1



moments=[ 0.0837  0.6737  0.0438  0.0698  4.4616];
 
actual_std=[ 0.0008  0.0042  0.0005  0.0010  0.1057]
actual_var= actual_std.* actual_std*1000;

 
m_target = moments';
m_scale=[actual_var'];

    weight_j=0.57;    



for jj=1:s;
    p0=param(jj,1);
    mu_c=param(jj,2);
    sig_eps_a=param(jj,3);
    rho_a=param(jj,4);   
    
    
   [f2, avfracup2,    med2,    interq2, kur2]=geNCalvoPlus_SMM(p0, mu_c, sig_eps_a, rho_a, weight_j);
 
   m_model=[f2; avfracup2;    med2;    interq2; kur2];
   G = m_target - m_model;
% Weighting matrix W
size_weights = max(size(m_target));
scale = (ones(size_weights,1).*[1;1;1;1;1])./m_scale; 
% take probability of default with 'double importance'
W = eye(size_weights).*(ones(size_weights,1)*scale');
clear scale size_weights

dist_sqrd = G.*G;
W_vector = [W(1,1);W(2,2);W(3,3);W(4,4);W(5,5)];
dist_weighted = dist_sqrd.*W_vector;
% Calculate new SMM estimator:

x = G'*W*G;

   
   
   
  mat_result19(jj,:)=[param(jj,:) x f2 avfracup2      med2    interq2 kur2  ]

  save mat_result19 mat_result19;
end


toc

disp('Program finished running on date/time:');
disp(datetime);
